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ABSTRACT 

We consider the problem of tidal disruption of stars in the centre of a 
galaxy containing a supermassive binary black hole with unequal masses. We 
assume that over the separation distance between the black holes the gravita- 
tional potential is dominated by the more massive primary black hole. Also, 
we assume that the number density of stars is concentric with the primary 
black hole and has a power-law cusp. We show that the bulk of stars with 
a small angular-momentum component normal to the black-hole binary orbit 
can reach a small value of total angular momentum through secular evolution 
in the gravitational field of the binary, and hence they can be tidally disrupted 
by the larger black hole. This effect is analogous to the so-called Kozai effect 
(Kozai 1962, Lidov 1961, 1962) well known in celestial mechanics. We develop 
an analytical theory for the secular evolution of the stellar orbits and calculate 
the rate of tidal disruption. We confront our analytical theory with a simple 
numerical model and find very good agreement. 

Our results show that for primary black- hole mass ~ 10 6 — 1O 7 M , the 
black- hole mass-ratio q > 10 -2 , cusp size ~ lpc, the tidal disruption rate 
can be as large as ~ 10~ 2 — lM Q /yr. This is at least 10 2 — 10 4 times larger 
than estimated for the case of a single supermassive black hole. The duration 
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of the phase of enhanced tidal disruption is determined by the dynamical- 
friction time scale, and it is rather short: ~ 10 5 yr. The dependence of the 
tidal disruption rate on the mass ratio, as well as on the size of the cusp, is 
also discussed. 

Key words: black hole physics - galaxies: nuclei 

1 INTRODUCTION 

The problem of the tidal disruption of the stars by a supermassive black hole has attracted 
the attention of many researches starting from the pioneering paper of Hills 1975, who 
proposed that the stellar gas released from the tidally disrupted star could power quasars 
and AGN's. Later, it has been realized that the tidal disruption event could manifest itself 
as a luminosity flare in the centre of a galaxy (e.g., Gurzadian & Ozernoi 1981, Lacy et al. 
1982, Rees 1988). The observational detection of these flares have recently been reported 
(e.g., Komossa et al. 2003 and references therein). 

In order to be tidally disrupted a star must have its orbital pericentre sufficiently close 
to the central black hole, and accordingly, its angular momentum must be sufficiently small. 
Thus, the presence of the central black hole give rise to a 'loss cone' in the distribution 
of stars in angular-momentum space (e.g., Frank & Rees 1976), where the number of the 
stars is quickly exhausted due to tidal disruption, and the process of refilling the loss cone 
determines the rate of tidal disruption. Estimates of the rate of tidal disruption are typically 
based on the assumption of two-body gravitational interactions between the stars as the 
dominant process of angular-momentum change. Since the relaxation time associated with 
this process is large for realistic models of galactic centres, these estimates give rather low 
values of the disruption rate, of the order of 10~ 7 to 10 _4 M o /yr depending on the properties 
of the stellar population of the galactic centres and the mass of the black hole residing in the 
centre (e.g., Frank & Rees 1976, Lightman & Shapiro 1977, Shapiro & Marchant 1978, Cohn 
& Kulsrud 1978, Magorrian & Tremaine 1999, Syer & Ulmer 1999, Wang & Merritt 2004 and 
references therein). It has been pointed out by Romberg 1968 and by Polnarev and Rees 1994 
that a secondary black hole of a smaller mass could intrude in the galactic centre as a result 
of the process of merging of two galaxies and perturb stellar trajectories, leading to non- 
conservation of the angular momentum and an enhancement of the tidal disruption rate (for a 
review of astrophysical processes occuring in galactic centres containing a supemassive black 
hole binary see e. g. Merritt & Milosavljevic 2004 and references therein). The purpose of 
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this paper is to calculate this enhanced tidal disruption rate in a simple model of interaction 
between the binary black hole and the stars. 

We assume that the binary black hole consists of two black holes of unequal masses 
moving around each other in circular orbits deep inside the potential well of the primary 
black hole. We take into account the gravity of the central stellar cluster in the form of small 
corrections to the motion of stars. It has been pointed out that when the gravity is dominated 
by the central black hole, the number density of stars increases toward the centre according 
to a power law and such a distribution of the stars has been called 'the cusp' (Peebles 
1972, Bachall & Wolf 1976, Young 1980). We use this power law distribution as the initial 
condition for our problem. We consider only stars with orbits within the orbit of the binary 
(which we will call the 'inner problem'). We show that in this case the main process leading 
to the replenishment of the loss cone is a secular effect increasing the eccentricities of the 
stars called 'the Kozai effect' (Kozai 1962, Lidov 1961, 1962). The qualitative explanation of 
the process is very simple. Let us assume as a zeroth approximation that the stars move on 
Keplerian orbits in the gravitational field of the primary black hole and that gravitational 
interaction with the secondary black hole introduces slow secular changes to the orbital 
elements. Neglecting the possibility of resonance interactions and considering the evolution 
of the orbital parameters on a time scale scale much larger than the orbital periods of the 
star and of the binary, it is well-known that the secondary black hole can be approximated 
by its orbital average, which is a mass ring. Since the gravitational field of a ring is static 
and axially symmetric, the energy and the angular-momentum component normal to the 
ring plane are conserved. However, in the course of the secular evolution the total angular 
momentum varies with time, and it can be as small as the constant normal component. Since 
the number of stars whose normal component of angular momentum lies within loss-cone 
size is much larger than the number of stars whose total angular momentum lies within the 
loss cone, the number of stars which can be tidally disrupted is significantly increased by 
the presence of the secondary black hole. 

In the next Section we formulate the model problem, which takes into account the secular 
evolution of the orbital elements of the stars from both the secondary black hole and the 
stellar cluster. The analytical solution of the secular problem is discussed in Section 3. The 
tidal disruption rate based on this solution is evaluated in Section 4. In Section 5 we present 
a numerical approach to the problem and check the validity of our analytical results. In 
Section 6 we use the results obtained in the previous Sections to estimate the enhanced tidal 
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disruption rate in astrophysical systems. This Section is written in a fairly self-contained 
way, so readers interested in astrophysical applications but not in the derivations of the 
relevant equations can go directly to this Section. Section 7 presents our conclusions and 
discussion. 

2 FORMULATION OF THE MODEL PROBLEM 

As explained in the Introduction, we would like to know how many stars are supplied on 
low angular-momentum orbits due to gravitational interaction with the supermassive binary 
black hole consisting of a primary black hole of mass M residing in the centre of a galaxy 
and an intruding secondary black hole with a lower mass m. The gravitational interaction of 
the binary black hole with stars of the stellar cluster is an extremely complicated physical 
process which involves dynamics of the stars in the field of two gravitating masses, effects 
of a slow change of orbital elements of the binary black hole due to dynamical friction and 
emission of gravitational radiation (e.g., Polnarev & Rees 1994, Quinlan & Hernquist 1997, 
Merritt & Milosavljevic 2004 and references therein) and due to interaction with an accretion 
disc (Ivanov, Papaloizou & Polnarev 1999), and effects caused by the gravitational potential 
of the stellar cluster itself. Obviously, this problem cannot be solved analytically. In order 
to make the problem analytically tractable while retaining the main features expected in an 
astrophysical situation, we consider a model where the orbital elements of a star are evolving 
in Newtonian gravitational field of two masses moving around each other on circular orbits. 
The separation distance D between two black holes remains constant with time. We also 
take into account the most significant effect of the gravitational field of the cluster on the 
stellar orbit, which is the apsidal precession of the pericentre of the orbit. The rate of this 
precession is calculated in Appendix A. Apart from the apsidal precession this model is the 
three-dimensional circular restricted three-body problem of celestial mechanics. 

Since we are interested in a statistical description of the stellar orbits, we assume below 
that the stellar orbits are characterised by phase space distribution function f(x,v,t). Now 
let us assume that at time t — the distribution function is not perturbed by the secondary 
black hole and equals its unperturbed stationary value fo(x,v). After this moment of time 
the gravitational influence of the secondary hole is turned on and the energies and angular 
momenta of the stars (defined in the frame centred at the primary black hole) start to evolve 
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due to perturbations induced by the secondary.* During the evolution the specific angular 
momentum of some stars gets smaller than a certain value It, corresponding to the tidal 
loss cone. These stars are assumed to be destroyed by the tidal forces of the primary black 
hole and are not considered further. 

Before going to details we would like to introduce some notation, discuss some additional 
simplifying assumptions used later on, and specify the initial conditions for the problem. 

2.1 Notation and additional simplifying assumptions 

When the semi-major axis of a star is sufficiently large, close encounters of the stars with 
the secondary black hole are possible. These encounters change the energy and the angular 
momentum of the stars in a very complicated way, and the analytical treatment of these 
encounters is rather difficult. In this paper we do not consider these complicated issues 
reserving them for future work and discuss only the 'inner problem' concentrating on the 
stars with semi-major axis a < D/2. In this case the orbit of the star does not cross the orbit 
of the secondary black hole, and the theory of evolution of the star's orbit is sufficiently simple 
(see below) . In the following discussion we introduce dimensionless variables associated with 
the binary orbit: a = a/D, L = Lj \J GMD, etc. In terms of these variables the problem 
depends only on two parameters: the dimensionless absolute value of the angular momentum 
associated with the tidal loss cone l T = L T /VGMD and the mass ratio q = m/M. Both 
parameters are assumed to be small, in the range q, It ~ 10~ 3 to 1CT 1 (see Section 6). 

2.2 The initial distribution function 

The number of stars supplied into the tidal loss cone obviously depends on the form of the 
initial distribution function Jq. In this paper we are interested in effects occurring deep inside 
the potential well of the primary black hole, where the number density of stars grows steeply 
near the centre of the stellar cluster, forming the so-called 'cusp' in the stellar distribution. 
In order to model this property of the central stellar clusters, we use the well known power- 
law isotropic distribution function of the stars in the phase space (e.g., Spitzer 1987, Bahcall 
and Wolf 1976, Lightman & Shapiro 1977) 

fo = CE p , 1 

* We neglect the change of the energies and angular momenta due to the usual gravitational scattering of the stars. This effect 
seems to be very slow comparing to the effects determined by the secondary black hole. 
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where C is a constant, E = GM/r — v 2 /2 is the binding energy of a star per unit mass in 
the field of the primary black hole of mass M, and v is the velocity of a star. According 
to theoretical estimates the parameter p < 3/2 should be rather small, and it can take 
two preferable values: p = (Peebles 1972, Young 1980) and p — 1/4 (Bahcall and Wolf 
1976). The value p = corresponds to adiabatic growth of the cusp in an unrelaxed stellar 
cluster containing a growing black hole and the value p — 1/4 corresponds to the situation 
when the characteristic relaxation time is smaller than the life time of the system. Since the 
relaxation time is typically rather large for the systems of interest we assume below p = for 
our estimates of the tidal disruption rate in astrophysical systems (see Section 6). However 
we present results from which the tidal disruption rate for the general case p ^ can be 
easily obtained (see Section 4). 

Also, let us note that after the binary black hole merging a significant part of stars 
may be expelled from the system and the cusp may be destroyed. In fact, the binary black 
hole scenario has been discussed to explaine rather shallow density profiles observed in the 
centers of giant central galaxies situated in the centeres of galaxy clusters (e.g. Faber et al, 
1997). We do not consider this more complicated situation in this Paper assuming that the 
binary black hole is formed the first time in a particular galaxy of interest. 

The number density of stars is obtained from (1) by integration over velocity space: 

r \ -(3/2+p) 



n[r) = n„ \ — 



where r a is the radius of the cusp and n a is the number density of stars at r a . The constant C 
in equation (1) is obtained by integration of equation (1) over velocity space and comparison 
the result of integration with equations (2). It is related to n a as 

n„ 

= 2 " {—) B(3/2; P+ 1)' 3 

where B(x,y) is a Beta function t The total mass of stars inside a sphere of radius r < r a 
is given as 

M (r) - 4tt mstUa r 3 ( T 4 

Mst{r) ~ An W^pj a w) ' 

where m s t is the star's mass and we assume hereafter that the cusp consists of stars of equal 
masses which are of the order of the solar mass M . The total mass of stars in the cusp must 
be of the order of the black hole mass M. In this paper we assume that these two masses 

t Note the useful relations: B(x, y) = B(y, x) = ^^f^ > wricre F(x) is the Gamma function. 
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are equal: M st (r a ) = M. This gives the normalisation condition: 

_ (3/2 -p)M 
4m st 7ir* 

The initial distribution of the stars over the binding energies can be obtained from equation 
(1) by integration over the surface in the the phase space corresponding to a given value of 
the energy E (e.g., Spitzer 1987): 

N (E) = C E E p ~ 5 / 2 , 6 

where Ce = V% (ttGM) 3 C. Below, instead of the distribution of the stars over the binding 
energy, it will be convenient to use the distribution over the semi-major axes a: 

N (a) = C a a l ' 2 - p , 7 

where C a = n 3 2 2 ~P {GM) 3/2+p C and we use E = GM/{2a). 

We also need the distribution of stars in inclination angle i and angular momentum. To 
obtain this function let us introduce the Cartesian coordinates with origin at the primary 
black hole and the z axis normal to the orbital plane of the binary black hole. The inclination 
angle i is measured from the z axis and the z component of angular momentum, L z , is related 
to its absolute as L z = Lcos(i). Assuming isotropy of the initial distribution function / and 
using equation (7) we immediately obtain the initial distribution over the semi-major axis, 
inclinations and the absolute values of angular momentum, No(a,i,L) 

iV (a,z,L)=iV (a)^^, 8 

where L a = \/GMa is the angular momentum corresponding to the circular orbit with 
radius a. The initial distribution function over semi-major axis, absolute values and the 
^-components of angular momentum follows from equation (8): 

N Q (a,L,L z ) = N (a) — . 9 
3 THE SECULAR EVOLUTION 

The bulk of stars supplied to the loss cone will have had small initial values of L z (see the 
next Section and Fig. 6 for numerical results.) This effect is determined by a particular 
feature of the secular evolution of the star's orbital parameters in the gravitational field of 
the binary. For the Keplerian inner problem it has been pointed out by Kozai 1962 and 
Lidov 1961, 1962 that although the a and L z are approximately conserved on time-scales 
much larger than the orbital periods of the star and the binary, while L varies with time. 
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Therefore, a star with small L z initially can, in the course of secular evolution, reach a 
small-enough value of L to be eventually tidally disrupted by the primary black hole. 

Kozai and Lidov independently derived the evolution equations for the orbital parameters 
using the methods of celestial mechanics. In general, these equations have a rather compli- 
cated form. However, they are simplified significantly when only the quadrupole contribution 
of the gravitational potential of the secondary is taken into account (e.g., Innanen, Zheng, 
Mikkola and Valtonen 1997, Kiseleva, Eggleton & Mikkola 1998). We use the evolution equa- 
tions derived in these papers and add an extra term describing the apsidal precession caused 
by the gravity of the stellar cluster. 

After doubly averaging over the periods of the star and the binary the relevant equations 
take the form 

Tk~t = -ee sin 2 (i) sin(2cj), 10 

OjV Zi 

di 5 

T K — = --e 2 sin(2i) sin(2cj)/e, 11 

Tr-j- = (2(1 - e 2 ) + 5 sin 2 (w)(e 2 - sin 2 (?)))/e - «e, 12 

Tjf = -(1 + e 2 (5 sin 2 H - 1)) cos(i)/e. 13 

UjL 

Here e, i, uj and \ are the eccentricity, inclination, the argument of pericentre and the 
longitude of ascending node, respectively, and e = \J\ — e 2 . T K gives the characteristic 
time-scale of the Kozai effect: 

3q V GM 

The last term in equation (12) describes the apsidal precession induced by the stellar 
cluster (see Appendix A). The parameter k is given by 

k = k ,^)(MM\=Ik(^1\. 15 



,P or b/ \ M J 3n \ m 
where P orb = 27ry / a 3 / (GM) is the orbital period of the star, M st (a) is the total mass of the 
stars inside a sphere of the radius a, and the positive coefficient K is calculated in Appendix 
A: 

In our estimates below we use K = 4a/2 corresponding to a cusp with the parameter p = 
(see equations (1-9)). It is very important to note that the apsidal precession caused by the 
stellar cluster has negative sign while the sign of the apsidal precession caused by the Kozai 
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effect is positive. This stems from the fact that the gravitational force acting from the cluster 
is always directed inward, to the centre of the cluster, while the averaged gravitational force 
acting from the secondary black hole in the inner problem is directed outward. The opposite 
directions of the gravitational forces result in the opposite signs of corresponding precession 
rates. 

The set of equations (10-13) has two integrals of motion: 

A = ecos(i), Q = e 2 (5sin 2 (i) sin 2 (cj) + re - 2). 17 

It is obvious that the integral A is proportional to L z . Using the expressions (17) we can 

obtain the evolutionary tracks of the dynamical system (10-13) on the plane (e,u>): 

6 2 ( g + (2 _ K)e 2 } 

8111 M = 5eV-A») ' 18 
Also, we can express the inclination i and the argument of pericentre uj in terms of the 

eccentricity e with help of the integrals (17), substitute the result in equation (12) and 

obtain a single evolution equation for the eccentricity e: 

T K ej t = ±^{((3 + k + Q - 5A 2 )e 2 - (3 + re)e 4 - Q){Q + (2 - re)e 2 )}. 19 

The solution of equation (19) can be expressed in terms of elliptical integrals (see below). 

It can be shown that in the course of evolution of the dynamical system (10-13) the change 
of the eccentricity with time is cyclic (we will call these 'Kozai cycles') and the argument 
of pericentre uj can either evolve monotonically or librate around the values uj = ±7r/2. 
Accordingly, solutions of equations (10-13) can be divided on two separate types of rotating 
(r) and librating (/) solutions. Let us discuss each types in turn. As follows from equation 
(10) the eccentricity e as a function of time can have its extremum only when uj = or 
7r/2. In the case of the solutions of (r) type both values of uj = and uj = tc/2 are possible. 
Let us call the values of eccentricity and inclination corresponding to uj = (uj = n/2) as 
eo and Iq (ei and The integrals (17) can be expressed either in terms of i$ and eo or %\ 
and e 1 . The corresponding expressions can be equated to each other and we can obtain the 
inclinations iq and i\ as functions of eo and e± from the resulting equations: 

sin 2 fa) = ^- ( 4 ~ ( 2 - « + (3 + «)e?), sin 2 (^) = (2 ~ ~ eg) . 20 



A 2 = ^^((3 + n)e\ + (2 - re)e 2 ), Q = -(2 - re)e 2 . 21 



The integrals A and Q can also be expressed in terms of e and e\. 

5e 2 

As follows from equation (20) when re < re cr i t = 2 the maximal (minimal) value of the 
eccentricity is attained at e = e\ (e = e ) and in the opposite case of a large re > re crit the 
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Figure 1. An example of Kozai cycles. The dependencies of eccentricity e (solid curve) and cos(i) (dashed curve) on time (in 
units Tfc). The minimal (maximal) eccentricity is equal to 0.5 (0.9999). The minimal value of cos(i) Ri (i ~ tt/2) corresponds 
to the minimal value of the eccentricity, and the maximal value of cos(i) 0.8 corresponds to the maximal eccentricity. The 
period of the cycles Px ~ l-5Tx is in agreement with its theoretical value (see equation (24)). The case shown corresponds to 
the Keplerian regime (k = 0) and the evolution of (r) type. 

maximal (minimal) value of eccentricity is equal to e (ei). In the last case the minimal 
value of the eccentricity t\ is always smaller than the critical value 



K 



^crit 



and the maximal value e is limited from above: e < min (l, yfrzfeij . It turns out that 
when k > K cr i t only a very small number of the solutions of the (r) type would lead to a 
significant increase of the eccentricity, and accordingly, to a significant decrease of the total 
angular momentum, during the secular evolution. Therefore, these solutions are not relevant 
to our problem and will not be considered further (see also the end of this Section). 

When k < K crit the solutions of the (r) type describe periodic changes of eccentricity and 
inclination, and the minimal (maximal) value of the eccentricity corresponds to the maximal 
(minimal) value of the inclination (see Fig. 1). 

Substituting expressions (21) in equation (19) and integrating the result we find the 
period of the Kozai cycles: 



Pk = 



^(2-«)(( 3 + «)e5 + (2-«)e§/e?) 
where K(fc) is the elliptic integral of first kind and 



K(y)T K , 



22 



(3 + «)(!- eg/eg) 
^ (2-«)eg + (3 + «)ef 



23 
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When k = and e <C e±, e 1 ~ 1 we have 

: / /a 4 \ 

Tx, 24 




where we use the asymptotic relation K.(k — > 1) ~ In y, and A;' = y/l — k 2 . It is interesting 
to note that when e ~ e\ 1 the period P° is finite: p» pa T^. From equation (22) 



it also follows that when k — > 2 the period increases as oc 1/y (2 — k). 

Now let us discuss the librating solutions where both minima and maxima of the eccen- 
tricity are attained at ou — ±7r/2. In this case the corresponding inclination angles and the 
integrals (17) can also be expressed in terms of the minimal and maximal eccentricities: 



1 

5 

and 



sin 2 (i ± ) = -(2 + 34 - «(1 ~4))' 25 



a2= (3 + «) (1 _ e 2 )(1 _ e 2 )j g =(3 + /t)e 2 e 2 ; 26 
o 

where, by definition, e_ (e+) is the minimal (maximal) value of the eccentricity and i± 
are the inclination angles corresponding to e±. From equation (25) it follows that when 



k > K crit = 2 only sufficiently large eccentricities are allowed: e± > e cr i t = y f^f- 

Similar to the previous case the expression for the period of the Kozai cycles is obtained 
by substitution of equation (26) in equation (19) and integration of the resulting differential 
equation. When k < k ci h we have: 



K 



e + y / (3 + /t)((3 + /t) e 2 + 2-k) 
where 



K(y)T K , 27 



(2 -«)(!- el/4) 28 
^ (3 + K)e 2 _ + 2 - k ' 

When k = 0, e + ~ 1 and e_ <C 1 we can obtain an approximate expression for very 
similar to the expression (24): 

F? ( 1% A \ 

29 




In the opposite case k > K CTit we have 

2 



K 



(3 + K)e-^/e 2 + -e 2 clit 
where 



K(y)T K , 30 



^crit 

y = 



1-1/4 31 



e- ^ 1 -4it/e 
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Figure 2. We show the evolutionary tracks of the dynamical system (10-12) on the plane (x\,X2) for the purely Kcplerian case 
(k = 0). The value of the integral A is taken to be 1 ■ 10 -3 for all curves. The solid curves represent the solutions of (r) type 
and the dashed and dot-dashed curves represent the solutions of (I) type. Different curves of the same type are characterised 
by different possible minimal values of the eccentricity, and respectively, by different maximal values of the 'radius' e. For the 
curves of (r) type the minimal value of the eccentricity eo is attained at ui = and ui = tt (when the curves cross the line 
X2 = 0), and we present the cases eo = 0.2, 0.5 and 0.8. The curves of (I) type have the minimal eccentricity e_ at u) = ±7r/2 
(when the curves cross the line xi = 0) and again we show the cases e_ = 0.2, 0.5 and 0.8. Two dotted lines separate the 
regions where the librating and rotating regimes of evolution are possible (see the text). 

Note that when k = n crit both expressions (27) and (29) give the same result: 

Pk = ■= T K . 

5e_e + 

When e_ — > e cv \ t the period diverges logarithmically. Assuming that 5 = (e_— e cri t)/e cri t 
1 we find from equation (30) 

^J5(k-2) V V (3 + ^)5/ 
where we also assume that e + 1 (see also equation 34). Note that the expression (32) is 
approximately valid only when 5 < <5* = ^ . Also, let us point out that when both e_ and 
e + are close to unity we can obtain another approximate expression: Pk ~ T Tk ■ 

v /5(k+3) 

3.1 The structure of the phase space 

It is convenient to represent the solutions of the different types on the plane (xi = ecos(cj), 
X2 = esin(u;)). Since e = \J(1 — e 2 ) < 1, all solutions are situated within the unit circle. 
Provided that the integral A is fixed the curves corresponding to different solutions belong 
to the same two-dimensional dynamical system and do not cross each other, and we use the 
small A = 1 ■ 10~ 3 hereafter. We show the solutions in Figs. 2-5. Fig. 2 and Fig. 3 represent 
the case of k < /t crit = 2, specifically k — and k = 1.9. The curves representing the 
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Figure 3. Similar to Fig. 2 but now k = 1.9 is just below the critical value re cr it = 2. All other parameters are the same as 
in Fig. 2. The shape of the curves of the different types is similar to the case re = 0, but the lines separating the regions of 
different evolution have smaller inclination with respect to the axis xi and the librating solutions occupy more space on the 
plane (xi,X2)- 

solutions of (r) type have a 'propeller' like shape, and the curves representing the solutions 
of (/) type have a 'bulb' like shape with 'cones' directed toward the origin of the coordinate 
system. When the eccentricity is not very large and the integral A is small, the curves 
separating the regions in the plane (xi,x 2 ), which correspond to the solutions of different 
types, are approximately two straight lines. Setting the value of eo to zero in equation (21), 
substituting the result in equation (18) and neglecting the value of A in denominator of the 
result we can find the inclination angles, o; C rit, of these lines 

sin (uait) = ±]j^-^. 33 

From equation (33) it follows that when k — > K crit , u crit — > 0. 

Fig. 4 and Fig. 5 represent the case of k > /t cr it- As we discussed above, in this case, the 
behaviour of the solutions of (r) type changes drastically. Contrary to the case k < /t crit , 
the solutions have minimal eccentricity at u — ±tt/2 and maximal eccentricity at uj = 0, 
7r. Also, the angle u> is rotating in the opposite negative direction compared to the previous 
case. The behaviour of solutions changes drastically and sharply when k exceeds /t crit . For 
example, in Fig. 4 we show the case of k — 2.1 which is just above the critical value of k — 2. 
One can see from this figure and Fig. 5 that only the curves of (r) type having the maximal 
values of eccentricity very close to the 'critical' value e CT \ t = \J can have a small value 
of the minimal eccentricity. Since we are interested only in the solutions which can have a 
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Figure 4. The case of re = 2.1 > re cr it is shown. The value of A is the same as in the previous cases. The shape of the curves 
of (r) type is drastically different from the cases re = and re = 1.9. Both rotating and librating curves have the minimal 
values of the eccentricity when crossing the line x\ = (a; = ±7r/2). When this value is smaller (larger) than the critical value 
e cr i t ~ 0.14 the curves belong to the family of (r) type curves ((I) type curves). As follows from the text for the curves of (r) 
type we define the value of eccentricity corresponding to u> = ±p/2 as e\. We show, respectively, the curves corresponding to 
ei = 0.05, 0.1, 0.12 (the solid curves) and the curves corresponding to e_ = 0.16, 0.2 and 0.5 (the dashed and dot dashed 
curves). The dotted 'critical' curve represents the boundary between the regions on the plane (xi,X2) where the motion of 
librating and rotating type is possible. 

very large value of the eccentricity the solutions of (r) type are not interesting for us when 
k > 2. On the other hand the librating solutions can produce a very large eccentricity, and 
therefore they are relevant to our problem even when k»1 (see Fig. 5). The solutions of 
both types are separated by a 'critical' curve which can be obtained from equations (18) 
and (26). For that we set e_ = e crit in equation (26) and substitute the result in equation 
(18). We have 

. . , U-2e 

sin(w crit ) = ±y — -■ 34 

From equation (34) it follows that when k increases the size of the regions on the plane 
(xi,X2) corresponding to the librating solutions decreases. Equation (32) also shows that 
the solution (34) formally has infinite period. 



4 THE SUPPLY RATE OF THE STARS TO THE LOSS CONE 

Now we would like to estimate the number of stars supplied to the loss cone as a function 
of time assuming that the initial distribution function is given by equation (9). To begin 
with we consider the Keplerian = 0. We also assume that the Kozai effect is solely 

responsible for the change of the angular momentum of the stars. As we have mentioned 
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Figure 5. Similar to Fig. 4 but corresponding to re = 10. In this case e cr i t as 0.784 and we show the curves corresponding 
to ei = 0.5, 0.75, 0.782 (the solid curves) and e_ = 0.82 and 0.9 (the dashed and dot dashed curves). Note that apart from 
a small region near the origin of coordinates the dotted curve representing the boundary between the corresponding regions 
almost coincide with the curve corresponding to ei = 0.782. 

above, the z-component of the angular momentum is approximately conserved in this case, 
and the number of stars per unit semi-major axis that can be supplied to the loss cone 
and tidally disrupted due to this effect can be estimated as the number of stars (per unit 
semi-major axis), Nq(o) with L z smaller than the size of the loss cone Lt- After integration 
over L z < Lt and L of the distribution (9) we obtain 

K( a ) = dL dL z N (a, L, L z ) = 2N (a)-^=, 35 
Jo J~l t va 

where N (a) is given by equation (7). 

Let us estimate the time dependence of the number of stars supplied to the loss cone. 

After the moment of time t — the number of stars with L z < Lt per unit of a, N c (a) is 

decreasing with time, and the rate of this decrease can be estimated as 

dN c (a) N c (a) 

where P K is given by equations (24) and (29) for the solutions of (r) and (1) type, respec- 
tively. As follows from this equation P K depends on the eccentricity e, and therefore, strictly 
speaking, in equation (36) we should divide the distribution function on P K before integra- 
tion of the distribution function over the absolute value of angular momentum. However this 
dependence is logarithmic, say P K (e = 0.5) ~ 1.487V and P^(e = 0.9) ~ Tk- Therefore, for 
simplicity, we neglect this dependence hereafter and set Pk = Tk from now on. Integrating 
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equation (36) we obtain 

N c (a) = iV c (a)e- (a3/2r) , 37 

where we use equation (14), r = ^-flbt, and Qb = \J^r- The number of stars iV*(a) supplied 
to the tidal loss cone and tidally disrupted follows from equations (35) and (37): 

N\a) = N c Q (a) - N c (a) = 2N (a)^=(l - e'^'^). 38 

Vfl 

The total number N l (t) of stars supplied to the loss cone is obtained by integration of 
equation (38) over a. The integration is straightforward, with the result 

N\t) = / o amaX daN\a) = j^C a D^a^{l - ^ll x -^h(2/3(l - p); *)), 

39 

where x = aj^f x r, and 7(/3, x) is the incomplete gamma function, a max is the upper limit of 
integration and we set the lower limit to zero. As we have discussed above we assume that 
5max ~ 1/2 for the inner problem^. Note, that to obtain equation (39) we use the explicit 
form of the initial distribution function over a given by equation (7). 

It is very important to point out that iV*(t) depends linearly on the small parameter 
l T - The fact that the dependence is linear leads to a significant increase of tidally disrupted 
stars in the presence of the secondary black hole (see Section 6). In contrast to this result the 
usual estimate of the tidal disruption rate due to two-body star-star gravitational scattering, 
in the so-called 'full' loss cone approximation (e.g., Hills 1975, Frank & Rees 1976) leads 
to a quadratic dependence. In the last case, the tidal disruption rate is proportional to the 
probability to find a star inside a circle in the angular momentum space of radius Lt, which 
is proportional to the area of this circle nL T . As we have mentioned above, in the presence 
of the secondary black hole only the z component of angular momentum is conserved, while 
total angular momentum varies due to the Kozai effect and can be evolved to a very small 
value L < L T . Hence, in our case the tidal disruption rate is proportional to probability of 
finding a star with L z < Lt- Note that in this respect our problem is similar to the problem 
of tidal disruption in a non-spherical galaxy where there is a tidal loss 'wedge' instead of 
the tidal loss cone, and the tidal disruption rate is enhanced (e.g., Magorrian & Tremaine 
1999). 

For comparison of our theory with results of numerical computations it is more convenient 
to use the probability Pr(t) for a star to find itself in the tidal loss cone. This can be obtained 

i Note, that the quadrupole approximation used in our analysis is, strictly speaking, invalid in the case of rather large a max = 
1/2. However, we assume that the estimates based on this assumption are correct in order of magnitude. 
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dividing the expression (39) by the total initial number of stars with semi-major axes smaller 
than a max , N tot . This number directly follows from equation (7) 

iVtot = —2^-al(t p . 40 



max 



(3/2 - p) 

Using equations (36) and (40) we obtain 

Prit) = { f- 2 ^ k (1 - fc^^-) 7 (2/3(l -p);x)). 41 

(1 -p) V a max 3 

The estimates shown above are approximately valid only for sufficiently large values of 
the loss cone It- The results of numerical calculations (see the next Section) show that when 
It is sufficiently small equations (39) and (41) significantly underestimate the number of 
stars supplied to the loss cone. To explain this effect let us consider the evolution of angular 
momentum on time scales of the order of the binary orbital period P5 = 27rf2^ 1 , which is much 
smaller than the Kozai time scale Tk- It turns out that the angular momentum oscillates on 
the time scale ~ P5 around an averaged value determined by the secular evolution on the 
time scale ~ T K . The largest amplitude of these oscillations is reached near the maximal 
value of the averaged eccentricity e max provided by the Kozai effect, and, respectively, near 
the minimal value of the angular momentum (see Appendix 1). The maximal amplitude, 
AL is estimated in Appendix B as 

A , 15 /3 ~ 2 

= ¥V l qa ' 42 

where Al = ALj\J GMD, and we take into account the fact that when k = both rotating 
and librating solutions have cos(i mm ) ~ (see equations (20) and (29)). It follows from 
equation (42) that when the size of the loss cone is very small, It <C q, the maximal value 
of L z of stars that can be supplied to the loss cone is of the order of Al 3> It- Accordingly, 

the number of stars supplied in the loss cone is determined by Al: 

r L a r AL a / 1^/3" 

A c (a) = / dL dL z N (a,L,L z ) = 2iV (a)-= = —J-qC a D 1/2 ~ p d 2 ~ p , 43 
Jo J-al \Ja 4 V 5 

where we use equations (7), (9) and (42). 

Similar to the case discussed above the number of stars per unit interval of a with 
L z < AL, N c (a) decreases with time. It is important to note that in the case of the small 
tidal loss cone the probability of a star with L z < AL to obtain the absolute value of 
angular momentum smaller than Lt during one Kozai cycle is of the order of Lt/AL 1. 
Accordingly, the rate of decrease of the number of stars with L z < AL (per unit of semi- 
major axis) can be estimated as 
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Integrating equation (44) we obtain 

N c (a) = N°(a)e-( Tl/VE \ 45 

where T\ = ^^Q^t, and the number of stars supplied to the loss cone is 

AT* (a) = iV c (a) - N a (a) = iV c (a)(l - e - (ri/v ^). 46 

The total number JV*(t) of stars supplied to the loss cone is obtained again by integration 
of equation (46) over a with the result 

N* = [ daN\a) = ^^^f^/^l -J), 47 

where 

and the probability for a star to find itself in the tidal loss cone, Pr x = N l /N tot is 



48 



Pr 15/3 (3-2p) 3/ 

m " 8 V5 (3-p) ^ maxU j - 



49 



The integral (48) can be estimated in terms of elementary functions in two limiting cases. 
Assuming that the argument in the exponent in (48) is sufficiently small, T\ < V^max, we 
can set a = a max in the exponent. In this case we have 

I w e-^i/v 7 ^^ 50 

and in the opposite case of a large time t\ > y/a^x, the integral is estimated to be 

7~2(3-p)^^e-( ri /^^). 51 

Taking into account the effect of additional oscillations of the angular momentum on the 
time scale of the binary orbital period significantly improves agreement between our theory 
and numerical calculations (see the next Section). 

4.1 The case k ^ 

When k 7^ the estimates of the number of stars supplied to the tidal loss cone must be 
modified. However, as we discussed in the previous Section, when k < n CT - lt = 2 the secular 
evolution of the orbital elements of the stars is very similar to the case k — 0, and therefore 
we do not discuss this case later on, assuming that our expressions (39), (41) and (47), (49) 
give a correct order of magnitude estimate. When k > 2 only the librating solutions of the 
secular equations (10-13) could produce a significant change of the eccentricity. Therefore, 
only the stars having orbital elements corresponding to the secular evolution of (1) type can 
be tidally disrupted provided that the L z of these stars is sufficiently small. As follows from 
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equation (34) the boundary between the solutions of (1) and (r) types can be defined as a 
closed 'critical' curve in the plane (e,u), and therefore, in order to estimate the number of 
stars going to the loss cone, we should introduce the initial distribution function of stars 
over the semi-major axes, the quantities e = y/l — e 2 , L z and the argument of pericentre 
u>, N (a, e, L z , uj). Taking into account that the distribution of stars is assumed to be the 
uniform over uj and that e = L/L a , the distribution function N (a,e, L z ,u) can be easily 
obtained from equation (9): 

N (a,e,L z ,u;) = ^-. 52 

We find the number of stars (per unit semi-major axis) that can be supplied to the tidal 
loss cone, Nq k , integrating the distribution function (52) over the surface under the 'critical' 
curve (34), and also integrating over L z : 

^o C ,» = f_[ dL z JJdS N (a, e, L z , u) = ^ arcsin (J^^j ^ N o( a )> 53 

where JJ dS (...) implies integration over the surface in the plane (e,u>) within the 'critical' 
curve (34), and U is the absolute value (in units of yjGMD) of the z component of angular 
momentum of the stars supplied to the loss cone. As we discussed above when It > g we 
have U ~ It and in the opposite limit 1* ~ A/ where A/ is given by equation (42). 

As follows from equations (30) and (32) the period of Kozai cycles Pk depends on k, 
and on the maximal and minimal eccentricities e±. Also, the period is formally infinite for 
the solution corresponding to the critical curve (34). We are interested in the case of a high 
maximal eccentricity e max pa 1 and the minimal eccentricity close to the critical eccentricity 
e cr i t . Therefore, we use the approximate value of Pk given by equation (32). Note that this 
value formally depends on the eccentricity e_ through the quantity 5 = (e_ — e cr i t )/e cr i t . 
However, since 5 enters the expression (32) only in the logarithmic factor we expect that 
the use of a definite value of 5 is irrelevant for our order of magnitude estimates, and for 



definiteness we set its value to the 'typical' value 5* = \J~-^- We have: 

2 In (2y/2) 2 
Pk ~ —f=====Tk ~ i T K . 54 

V 5 ( K - 2 ) V^-2) 

It is very important to note that both the number of stars supplied to the loss cone N£ (a) 
and the period of Kozai cycles P K are proportional to the same quantities corresponding to 
the case k — 0: Pk = C±Tk, and Nq k (o) = ^iV^a), where C\ is determined by equation 
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(54), and C 2 follows from equations (35), (43) and (53): 

C 2 = — arcsin ( J — - — | . 55 
7T \ V k + 3 

That means that the results obtained above for the case n = can be used when calculating 
the number of stars supplied to the loss cone for a non-zero value of k. We should only 
multiply the expressions (39), (40) and (47), (49) by the correction factor C 2 and also 
make an appropriate change of the time variables entering these expression with the help of 
equation (54). When l T ^> q we obtain from equations (41), (54) and (55) the expression for 
the probability for a star to find itself in the loss cone, Pr K (t): 



56 



where 

x K = x/d = - ^5(K-2)~a 3 ^ x qn b t. 57 
In the opposite case It <C q we use equations (49), (54) and (55) 



PnAt) = C 2 P ri (t) = ^ arcsin f ) ^(1 - 7(r 1>(C )), 58 



ti,« = y gTi/Ci = J ^r—hQbt, 59 



where 



and the quantity J is given by equation (48). Note that we discard the factor y | in equation 
(42) and accordingly in equations (44)-(49) to obtain equations (58) and (59). This accounts 
for the coefficient in equation (59). As follows from equation (25), for the eccentricities 
e_ ss e crit , sin(? min ) s^ and accordingly cos(? min ) s^ 1. This is different from the case of 
k = where cos(i min ) s^ y^. 

4.2 Characteristic rates and time scales 

In general the expressions for the number of stars supplied to the loss cone given by equations 
(39), (41), (47), (49) and (56), (58), are rather complicated functions of time. They also 
involve many parameters determined by the properties of a particular stellar cluster and a 
binary black hole. These complicated expressions are rather difficult to use. Also, in fact, 
they are not needed for rough estimates of the tidal disruptions rate in the presence of the 
black hole binary. Here, we would like to obtain much simpler order of magnitude expressions 
which will be used in Section 6 for these estimates. For simplicity, we also assume that the 
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mass of the stellar cluster is much larger than the mass of the secondary black hole, and 
accordingly, formally consider the case of high values of k ^> 1 . 

In order to obtain the total probabilities for a star to be tidally disrupted let us consider 
the limit t — > oo in equations (56) and (58). In this limit, the second terms in the brackets 
in (56) and (58) are equal to zero and we have 

x 4^(3-2^) l T 
Prtot = Pr K (t - oo « H~^^= 60 

7T (1 -p) VflmaxK 

from equation (56), and 

n / 15 v 7 ^ (3 - 2p) a 3 / 2 

Pr M = Pr lA t - oo) « _«^_W^ 9 , 61 

from equation (58). The characteristic tidal disruption rate Pr per star can be obtained 
from equations (56) and (58) by differentiation these equations over the time and setting 
t = in the resulting expressions. We have 

Pr = Pr K (t = 0) = Pr liK (t = 0) = ^ ^g o^gi^ 1 . 62 

Let us note that the quantity Pr does not depend on The characteristic time scales 
associated with the process of tidal disruption, T* and T* i can be obtained from equations 
(60), (61) and (62). When l T > q we have 



and in the opposite case I? < g we have 



V5 (5 - 2p) aj, 



«V2 



T„ 1 = F rM , 1 /Pr=Y^-^^ f (W. 

As we discussed in Section 2, the typical values of the parameter p are close to zero. Therefore, 
we can further simplify expressions (60)-(63) setting p = and a max = 1/2 and expressing 
k in terms of the ratio of the total mass of the stars to the mass of the secondary, M st /m 
with help of equations (15) and (16): 

SV2M st (a max ) 4 M st (D) M st (D) 

k = = — - ps 0.4 — -, 65 

3n m 3n m m 

where we take into account that when p = 0, M st (a) oc a 3 / 2 (see equation 4). 
We have 



9 

Pr = —ql T Q- 1 ^ \Aql T Vt b 66 

2"7T 



from equation (62), 



Prtot = J™ ™i T ~ l8.5 A /^-ir, 67 
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from equations (60), (63), and 

„ 15/15 / m nn I m 

from equations (61) and (64). Equations (66-70) are used in Section 6 for estimates of 
the tidal disruption rate in astrophysical systems. Although equations (66-70) have been 
formally obtained in the limit k ^> 1 they can also be used for order of magnitude estimates 
even when k ~ 1. 



5 NUMERICAL CALCULATIONS 

Now we would like to discuss our numerical calculations of the supply rate and compare them 
to the analytical theory developed above. In our numerical work we neglect the influence of 
the stellar cluster on the orbits of the stars, and accordingly set k = reserving the general 
case k 7^ for a future investigation. We consider evolution of a large number of particles 
in the Keplerian gravitational field of the binary for a sufficiently long time, assuming that 
all particles which have the absolute value of angular momentum smaller than a certain 
value l T are 'tidally disrupted' by the primary black hole. The 'tidally disrupted' particles 
are removed from the calculations. The initial orbital elements of the particles are randomly 
distributed according to the distribution (8) with p = 0. The initial semi-major axes of the 
particles are in the range (0.1 < a < 0.5). The calculations depend on two dimensionless 
parameters: the mass ratio q and the dimensionless size of the loss cone It, and we consider 
the cases q=l- 10^ 2 , 3 • 10~ 2 , 1 • 10" 1 and l T = 1 ■ 10" 3 , 3 • 10" 3 , 1 • 10" 2 , 3 • 10~ 2 , 1 • 10" 1 , 
respectively. For given values of q and It we evolve approximately 10 4 — 5 • 10 4 particles with 
different orbital parameters during the time t < t max = — fi^ 1 . 

In Fig. 6 we show the distribution of the 'tidally disrupted stars' over the initial values of 
the z component of angular momentum, for the case q — 0.1 and different values of It- One 
can see that that this distribution is strongly peaked at L z = 0. This strongly suggests that 
the effects of the Kozai type are responsible for the supply of the particles on low angular 
momentum orbits in our problem. The spread of the distribution decreases with decrease 
of It- However, even for the smallest value of It — 1 • 10 -3 the spread is of the order of 
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Figure 6. The distribution of the 'tidally disrupted particles' over the initial values of the z component of angular momentum 
(in unit y/GMa). The case q = 1 is shown. The solid, dashed, dotted, dot-dashed and dot-dot-dashed curves correspond to 
l T = 1 . l(r 3 , 3 • 1CT 3 , 1 • 10" 2 , 3 • ltr 2 , 1 • 10" 1 , respectively. 




q£i„t 



Figure 7. The dependence of the number of stars supplied to the loss cone (in units of the total number of stars participating 
in the calculation) as a function of time (in units of q~ 1 Q~^ 1 ). The size of the loss cone is large: It = 0.1. The three solid curves 
correspond to q = 1 • 10~ 2 , 3 • 10~ 2 , 1 • 10 -1 , with smaller values of ordinates (at a given moment of time) corresponding to 
smaller values of q. The dashed curve represents the theoretical solution (41) which docs not depend on q in these time units. 
The dotted curves represent the sum of the solutions (41) and (49). The smaller ordinates correspond again to the smaller 
values of q. 



~ 0.1 ^> It- This can be explained by the additional oscillations of the angular momentum 
due to gravitational influence of the secondary on the time scale of order of the binary orbital 
period, see equation (42) and Appendix B. The distributions corresponding to the smaller 
values of q have a similar character and are not shown here. 

In Fig. 7 we show the numerical and theoretical values of probability for a star to find 
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itself in the loss cone as a function of time, for a large value of the size of the tidal loss cone. 
As follows from the theoretical solution (41), in the simplest approximation this probability 
depends on the mass ratio only in the combination qVL^t. That means that when the size of 
the loss cone is large and fixed, the numerical values of the probability (solid curves) should 
be approximately the same being plotted as functions of the time variable qVL^t. In order to 
check this theoretical prediction we show the results of the numerical calculations for a very 
large and somewhat unrealistic size of the loss cone It = 0.1. As follows from Fig. 7 this 
theoretical prediction is qualitatively confirmed in the calculations, although the difference 
between the numerical curve corresponding to q = 0.1 and q = 3 • 10~ 2 and 1 • 10~ 2 is rather 
large, of the order of ~ 25-30 per cent. The agreement between the theory and numerical 
calculations can be improved significantly if we take into account the effect of oscillation of 
the angular momentum on the orbital time scale and add the probability Pr\ corresponding 
to this effect (see equation 49) to the probability Pr determined by the Kozai effect (equation 
25). § The curves representing the sum Pr tot = Pr + Pt\ are shown as the dotted curves in 
Fig. 7. The curves corresponding to the mass ratios q = 1 • 10~ 2 and 3 • 10~ 2 are close to 
the curve representing Pr(t). The difference between these curves and those which represent 
the numerical calculations is about ~ 10-15 per cent. The dotted curve corresponding to 
q = 0.1 is rather close to the numerical curve with the difference of the order of ~ 6-10 per 
cent. Therefore, when the size of the loss cone is large the use of the total probability Pr to t 
instead of Pr improves the agreement between the theory and the numerical calculations. 

In Figs. 8, 9 and 10 we show the dependencies of numerical and theoretical values of the 
probability for the particle to be supplied to the loss cone calculated for realistic values of 
the loss cone l T — 1 • 10~ 3 and l T — 3 • 10~ 3 . The mass ratio is equal to 0.1, 3 ■ 10~ 2 , 1 ■ 10~ 2 
in Fig. 8, Fig. 9 and Fig 10, respectively. The general behaviour of the curves is similar in all 
cases. The curves corresponding to It — 3 • 10~ 2 are better described by the theoretical curve 
representing the total probability Pr tot (t). On the other hand the curves corresponding to a 
very small It — 1 • 10~ 2 are much better described by the curve representing the probability 
Pri{t). Note, that the agreement between the theoretical and numerical values is very good 
in the case q — 0.1, and the theoretical and numerical curves almost coincide in this case. 
The disagreement grows with decrease of the mass ratio q, and in the case q — 1 • 10~ 2 is 

§ Note that we use the approximate values of the integral I in the expression for the probability Pri. When the time variable 
n < 6 wc use the expression (33), and when ti > 6 we use the expression (34). When ri = 6 both expressions are equal. 
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Figure 8. The dependence of the number of stars supplied to the loss cone (in units of the total number of stars) as a function 
of time (in units f^ 1 ) for q = 0.1 and the small values of the loss cone. The upper solid curve corresponds to It = 3 • 10~ 3 
and the lower solid curve corresponds to It = 1 ■ 10~ 3 . The dashed and dotted curves represent the theoretical values of Pri(i) 
and Pr tot (t), respectively. The upper (lower) dashed and dotted curves correspond to l T = 3 • 1CT 3 (It — 1 ■ 10~ 3 ). 
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Figure 9. Same as Fig. 4 but calculated for q = 3 • 10 2 . 

about ~ 10 per cent. Thus, in the case of a very low size of the loss cone the probability for 
the particle to be supplied to the loss cone is mainly determined by the effect of oscillations 
of the angular momentum on the binary orbital time scale, and respectively, by equation 
(49). In fact, the numerical results demonstrate that for a rough estimate in the case of 
a large ratio It/q we can use the expression for the probability determined by the Kozai 
effect (equation 41) while in the opposite case of the small ratio lx/q the expression (49) 
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Figure 10. Same as Fig. 4 and Fig. 5 but calculated for q = 1 • 10 



determined by the effect of the oscillations of the angular momentum. As we have discussed 
above, this conclusion is well motivated from theoretical point of view. 



6 ASTROPHYSICAL APPLICATIONS 
6.1 The size of the tidal loss cone 

We do not concentrate in this paper on complicated issues related to the problem of tidal 
disruption of a star, and use the simplest criterion of the tidal disruption event assuming 
that all stars having pericentres R p smaller than the so-caller 'tidal' radius Rt are tidally 
disrupted: 

) R st « 1.5 • 10 13 M 7 1/3 cm, 71 

m st J 

where we assume that all stars have solar masses and radii R st ss 7 • 10 10 cm and M 7 = 
M/1O 7 M . A simple estimate (Hills, 1975, Beloborodov et al 1992) shows that when the 
mass of the central black hole is larger than the 'critical mass' M ~ 1 • 1O 8 M the stars are 
swallowed by the black hole without disruption, and therefore the range of the black hole 
masses interesting for the tidal problem is limited from above: M < MO 8 M .H O n the other 
hand, at the present time there is no observational evidences of black holes masses within 

^ Note, however that a rotating black hole of mass ~ 10 9 Mq can actually disrupt the stars due to effects of General Relativity 
(Beloborodov ct al 1992, Ivanov, Chcrnyakova & Novikov, 2003). The tidal disruption rate is suppressed in this case by a 
geometrical factor, and wc neglect this effect later on. 
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the intermediate mass range M < 1O 5 M (e.g., Kormendy & Gebhardt 2001 and references 
therein), and therefore we assume that our parameter q cannot be smaller than ~ 10~ 3 . The 
specific angular momentum Lt associated with Rt is obtained from the standard Keplerian 
relation L T = ^2GMR T . Following from the previous Sections, in our calculations we use 
the dimensionless angular momentum variable It = L T /y/GMD : where D is the separation 
distance between the two black holes. The numerical value of It can easily be obtained from 
equation (71): 

l T « 3 • 10- 3 M 7 1/6 (/3ri) _1/2 , 72 

where D = D/r a and r\ = r a /lpc. We see a rather weak dependence of the numerical value 
of It on the parameters, especially on the mass of the central black hole. 

6.2 The dynamical friction 

When the secondary black hole enters the cusp it interacts with the stars in the cusp and 
continues to spiral down to the centre of mass due to the effect of dynamical friction. As 
a result the separation distance between the primary and the secondary black holes is a 
decreasing function of time D(t). To estimate the characteristic time scale associated with 
this process we have to specify the distribution function of the stars in the cusp as well 
as the orbital parameters of the binary. As we have mentioned above we assume that the 
binary is circular and use the isotropic distribution function (1) with p — 0. In this case the 
evolution law of the binary separation distance is especially simple: 

D(t)=r a exp(-t/T d . f .), 73 

where r a is the radius of the cusp (see Section 2), and the characteristic evolution time scale 
Td.f. can be easily obtained from the expressions given in Chandrasekhar 1960: 

Td.f. = ^K'Pa - 2 • 10 4 {^f) M 7 - 1/2 rf yr, 74 

where q = m/M is the binary mass ratio, A c ~ 15 is the so-called Coulomb logarithmic 
factor, r\ = r a /lpc, and P a is a typical period of the binary with separation distance D ~ r a 

P a = » 3 • 10 4 Mf 1/2 rf 2 yr. 75 

The evolution law (73) tells that the characteristic evolution time scale D/D is independent 
of the distance D and is always approximately equal to T d f. 

Equation (73) has its limits of applicability. When the binary evolves to the radius 
r m = r a q 2 / 3 where the total mass of the stars inside the binary orbit is of order of the 
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secondary mass m the evolution is halted (Begelman, Blanford & Rees 1980) Jl So we have 
D > g 2 / 3 . On the other hand, the expression (73) is approximately valid only when the time 
Td.f. is larger than the binary orbital period 

P b = 27Tfi 6 - 1 = P a D^\ 

hence D < ( 3 ^ A . When the mass ratio of the binary is larger than a certain value 
g crit both conditions are incompatible with each other and our simple description of the 
dynamical friction becomes qualitatively invalid. The value q CTl t can be found by equating 
Tdj. and Pb at the radius r m 

where we recall that M st = MD 3 ^ 2 (equations (4) and (5) with p = 0) and use A c = 15. 
When the mass ratio q is in the range: ~ 1 ■ 1CT 2 < q < q crit ~ 0.1 equation (73) is 
approximately valid only for the dimensionless separation distance D smaller than 

, 4/3 / \2/3 

(3rrA c g) 2/3 Q 2/3 \ Q J 

6.3 The tidal disruption rate 

The expression for the tidal disruption rate, M t is determined by the ratio of the character- 
istic time scale of the secular evolution given by equations (68) and (70) to the time scale 
of the dynamical friction Td.f. When this ratio is small, for a given value of the separation 
distance D almost all stars with semi-major axes of the order of D and sufficiently small 
values of the z component of angular momentum have a chance to be tidally disrupted. In 
this case M t can be expressed in terms of the time scale Td.f, the probabilities Pr tot and 
Pr totj i given by equations (67) and (69), and the total mass of the stars within the sphere 
of radius D, M st (D). When l T 3> q we have 

M t(D) „ «5**-C2, 

J-d.f. 

and in the opposite case ^ C 5 we should use Pr tot x instead of Pr tot in equation (78). When 
the time scales of the secular evolution are larger than T d j. the bulk of the stars have no 
time to be tidally disrupted, and the tidal disruption rate is determined by equation (66): 

M t (D) w PrM st (D). 79 

II In fact, the evolution may proceed due to the effect of gas dynamic drag exerted on the binary by e.g., an accretion disc 
(Ivanov, Polnarev & Papaloizou 1999). We do not consider this possibility in the present paper. 
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As follows from equation (72) the size of the loss cone increases when the separation 
distance D decreases. That means that, in principal, the ratio R = It/q is a function of 
scale. However, it is easy to see from equation (72) that when q > 1 • 10~ 2 we have R < 1 for 
the all scales of interest r m < D < r a , and when q < 3 ■ 1CT 3 we have R > 1 for the relevant 
scales. For simplicity we are going to consider separately two possible cases: 1) the case of 
a large mass ratio q > 1 • 1CT 2 and 2) the case of a small mass ratio q < 3 • 1CT 3 . 

In the case of the large mass ratio it can be easily shown that the time scale 7^./ is always 
smaller than or of the order of the time scale T^i given by equation (70). Accordingly, we 
can use equation (79) and obtain the explicit expression for M t with help of equations (66), 
(72) and (75): 

M t « 4 • 10- 1 ( V 1/3 M 7 5 /Vr 2 ^, 80 
\10~V yr 

where we use the radial variable ji = M st (D)/m instead of D. It is obvious that the variable 
ji takes its values in the range 1 < (j, < q" 1 . M t decreases with increase of \i and therefore 
the innermost part of the cusp /i = 1 gives the highest tidal disruption rate. Note, however, 
that this dependence is rather slow M t oc /i -1 / 3 . 

In the case of the small mass ratio there is a scale /z* where the secular time T* given by 
equation (68) is equal to T d f. 

/3-10- 3 \ 2 

/i* w I I . 81 

The tidal disruption rate is given by equation (80) for //>//*. In the opposite case we use 
equation (78) with the result 

Mt « 2 ■ 10- ( j^) 5/3 /! 1/6 M 7 5/3 r7 2 ^. 82 



yr 



The tidal disruption rate is maximal at \x = //*: 

A/, m «2.10->(^) 4/ VVr^f. 83 
Estimates of the tidal disruption rate in the astrophysical systems containing a single 
supermassive black hole give the value of the rate smaller than or of the order of M* ~ 
10 _4 M o /yr (e.g., Magorrian & Tremaine 1999, Syer & Ulmer 1999 and references therein). 
Equating M* to M max we can find the minimal mass ratio, g m j n sufficient for the tidal 
disruption rate determined by the presence of the binary to surpass the standard value 

q > q min » 10- 5 r 3/2 Mf 5/4 . 84 

There is another important limitation on our formalism determined by the relativistic 
effect of Einstein apsidal precession of the stellar orbits (Hopman & Alexander 2005). To 
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estimate the importance of this effect we use the standard expression for the change of the 
argument of peri centre due to this effect: 

6nGM 

0Ul re l — - » , 

c 2 ej,a 

where = It/V& is the value of the parameter e associated with the 'tidal' angular mo- 
mentum It- We calculate the corresponding precession rate (o re \ = 5u re i/P or b, and compare 
it with the precession rate determined by the Kozai effect in the limit of high eccentricity: 
io K ~ 5/(T K e T ) (see equation 12). We find that the Einstein precession is not important 
provided that the mass ratio is sufficiently high: 

M 5/6 

q > q re i ~ 1.5 • 10~ 3 ^=, 85 

where we set a = D/2 for a 'typical' stellar orbit.' 

Finally, let us estimate the duration of the phase of the enhanced tidal disruption rate 
due to the presence of the secondary. The corresponding time scale, T T follows from equation 
(73) 

T T = In (r a /r m )T dJ . = ^ In (g" 1 )^./. 86 

We also neglect the difference between D crit (equation 77) and r a which enters in the ex- 
pression under the logarithm. Substituting equation (74) in (86) we find that TV is rather 
short: 

T T = 6 • 10 4 ( 1 + -!- In — ^ Mf 1/2 rf V- 87 
V 4.6 q J 

7 CONCLUSIONS AND DISCUSSION 

We have shown in this paper that the tidal disruption rate in a galactic cusp containing a 
circular black hole binary can be as large as ~ O.l-lM /yr for a cusp size of the order of lpc. 
Cusps of this size have been observed in e.g. our own Galaxy and M32 (e.g. Kormendy & 
Gebhardt 2001 and references therein). In principle, the rate of gas supply to the centre of the 
galaxy is sufficient to power a quasar or a powerful AGN. However, the estimated duration 
of the stage of the enhanced tidal disruption rate ~ 10 5 yr seems to be much shorter than the 
expected life time of the central engine in quasars. Also, note that the tidal disruption rate 
sharply decrease with increase of the cusp size r a : oc r\ 2 , where r\ = r a /lpc. The obsevations 
suggest that the cusp size grows with mass of the central black hole with the typical sizes 
in the range 10 — 30pc for the masses of central black holes of the order of 10 7 — 1O 8 M 
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(e. g. Kormendy & Gebhardt 2001). So for an extended cusp of size ~ lOpc we get rather 
moderate values of the tidal disruption rate of the order of 1 ■ 10~ 3 — 1 • lO~ 2 M /yr. 

The tidal disruption events in a system containing a supermassive black binary black 
hole could have observational consequences on its own. A single tidal disruption event could 
manifest itself as sharp rise of luminosity of a galactic centre on time scale ~ lyr (Gurzadian 
& Ozernoi 1981, Lacy et al 1982, Rees 1988, Komossa et al 2003 and references therein). 
The repeating luminosity flares in a particular galaxy could suggest that there is a binary 
black hole in the centre. On the other hand, one can expect the increased probability of the 
tidal disruption events in the supermassive binary black hole candidates (e.g., Sillanpaa et 
al 1988, Valtaoja et al 2000). 

We assume that the stellar distribution has a pronounced density increase toward the 
centre. This assumption is well motivated from theoretical point of view and also has been 
confirmed by analysis of the distribution of stars inside our own Galaxy (e.g., Alexander 
1999) and some other galaxies (e.g., Faber et al. 1997). However, after the intervention of 
the secondary black hole and coalescence with the primary, a large fraction of the stars is 
expelled from the centre, and the density increase is smeared out. This may explain rather 
shallow density profiles observed in the centers of bright galaxies (e.g., Faber et al. 1997) 
Thus, subsequent coalescence events of two black holes within the centre of the same galaxy 
would be accompanied by smaller numbers of tidally disrupted stars. 

We have considered only the 'inner problem' where all stars going to the loss cone were 
assumed to be inside the binary orbit and showed that for the inner problem the bulk of 
the stars going to loss cone would have had a small value of the projection of the angular 
momentum on the axis perpendicular to the binary orbital plane before intervention of 
the secondary in the cusp. The mechanism responsible for the supply of the stars into the 
loss cone is analogous to the well known Kozai effect in celestial mechanics. However, the 
important effect of the oscillation of the angular momentum on the time scale of the order 
of the binary orbital period should also be taken into account in order to obtain the correct 
expression for the tidal disruption rate in a binary with a sufficiently large mass ratio. Also, 
we take into account the apsidal precession determined by gravitational potential of the 
stellar cluster itself and develop the analytical theory of secular evolution in the combined 
gravitational field of the binary and the stellar cluster. This theory might be useful in other 
astrophysical applications, such as e.g., the theory of dynamics of the triple stars (e.g., 
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Kiseleva, Eggleton & Mikkola 1998 and references therein), and the theory of dynamics of 
extra-solar planets around the binary or a single star (e.g., Papaloizou & Terquem 2001). 

The tidal disruption rate may be significantly enhanced if we take into account 'the 
outer' stars with semi-major axes larger than the half of the binary separation distance. 
The evolution of these stars in the binary gravitational field is mainly governed by two 
processes: the secular non-resonant interaction with the binary similar to the interaction 
discussed in our paper, and the effect of close encounters with the secondary which change 
the orbital elements of the star in a complicated chaotic way. Unfortunately, this problem 
looks very complicated for a thorough theoretical treatment and numerical calculations are 
very important in this case. 

Our estimate of the rate of supply of the stars to the tidal loss cone relied upon the 
effect of the secular increase of the eccentricity. We have analysed this effect in an idealised 
situation where only the Keplerian gravitational field of the binary and corrections due to 
the stellar cluster have been taken into account and the binary has had a circular orbit with 
the fixed separation distance between the primary and the secondary. In a real astrophysical 
situation the binary could have a substantial eccentricity (e.g., Quinlan & Hernquist 1997, 
Aarseth 2003) and the orbital parameters of the binary are evolving with time due to the 
effect of dynamical friction. Therefore, a very important question arises: to what extent is 
the effect of the secular increase of the eccentricity of the stellar orbits sensitive to possible 
corrections due to non-zero eccentricity and to a change of the orbital elements of the 
binary? First we would like to point out that when the eccentricity of the binary is not very 
large the effects determined by the eccentricity seem to be rather insignificant for the inner 
problem considered above. Indeed, in the quadrupole approximation used in our analysis the 
gravitational field of the eccentric binary is formally equivalent to the gravitational field of 
the circular binary provided that we use the parameter of the binary orbit ps = «b(1 — e%) 
instead of the separation distance D in equation (14) defining the characteristic time-scale 
of the Kozai effect, Tk (e.g Innanen, Zheng, Mikkola and Valtonen 1997). Of course, in the 
case of a highly eccentric binary the influence of the next order corrections in the small 
parameter a/D would be more prominent, but we expect that the qualitative character of 
the secular evolution would remain unchanged until the condition a/D < 1 holds for the 
stellar orbits of interest. When the quadrupole approximation is valid, a slow change of 
the orbital parameters of the binary also seems to be unimportant. In this case this time 
dependence enters the equations of the secular evolution (10-13) only through the time 
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dependence of the quantity T^(t) in front of the time derivatives. Introducing the new time 
variable t± — J dt/Tx we can eliminate this time dependence, and the secular evolution is 
the same in terms of the new time variable t\ as it was in the case of the non-evolving binary 
orbital elements. 

In this Paper we consider the black holes as point like sources of Newtonian gravity. 
However, as we discussed in the previous Section, the relativistic effect of Einstein apsidal 
precession (and, possibly, the effect of gravitomagnetic precession) could be important. A 
study of these effect could be a very interesting problem deserving a separate study. 

In conclusion, we would like to point out that for a quantitative calculation of the tidal 
disruption rate our approach is perhaps too oversimplified, and a thorough numerical study 
is needed. However, our model can be used for the order of magnitude estimates discussed 
above. It can also be used for guidance on the choice of parameters in numerical studies. 
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APPENDIX A: THE PRECESSION RATE 

Here we derive expressions for the apsidal precession rate of a stellar orbit induced by the 
gravitational field of the stellar cusp. As follows from the discussion above the mass of the 
stars in the cusp M st and accordingly the gravitational potential <& st depend on the radius 
r as a power of the radius 

/ r x 3/2-p Q M / r x 1/2-p 

M st — M ( — J , = (L\ , Al 

\r a J (1/2 -p)r a \r a J 

Assuming that r C r„ we can consider $ st as a perturbation. In this case the change 5oj of 
the angle ui per one orbital period is given by the expression (e.g., Landau, Lifshitz 1969) 

5u = 24 7 (t f dcf>(r 2 <S> st )) , A2 



dL \L Jo 

where the radius r is determined by the standard Keplerian relation 

r — A'i 

GM(l + ecos(0))' 

Substituting (Al) in (A2) and introducing the new variable e = \J (1 — e 2 ) instead of the 
eccentricity e we obtain 
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where s = 3/2 — p, a is the semi-major axis, we use and 



J= / rrr, AA 

Jo (l + ecos0) s+1 

When s ^ n/2 where n is an integer, the integral J can be expressed in terms of the 
Legendre function P s (x) as 

J = Tre-^P^e- 1 ), A5 
and P s (x) is expressed itself in terms of the hyper-geometric function F(a; /3; 7; x) as 
P.(x) = V ^^^ x'Fdl - s)/2; -s/2; 1/2 - S ; 

I^T^T372T^*" F(s/2 + 1; (s + 1)/2; s + 3/2; 1/l2) ' M 

where T(x) is the gamma function. 

We are interested in s = 3/2 corresponding to p — 0. In this case the integral J can be 
expressed in terms of the complete elliptic integrals as 

When eccentricity is small, the expressions are simplified significantly In the leading 
order in the small parameter e we obtain from equations (A3), (A5) and (A6) 

5u = -K (-Y e, A8 



where 

k — .H, 

(2s -1) T(s) 

Substituting p = 3/2 — sin equation (A9) we obtain 
In the case p = (s = 3/2) we use equations (A3) and (A7) to obtain 

/„\3/2 

5lo = -4V2 -) e. All 



k = ^,^,r( 8 _ + i/2)^ M 



A10 



r, 



Note that equation (All) follows from equations (A8) and (A9) when s = 3/2. 



APPENDIX B: EVOLUTION OF THE ANGULAR MOMENTUM 

In the quadrupole approximation the rate of change of square of the angular momentum of 
a star on an elliptic orbit due to gravitational perturbation of a perturber of mass m is given 
by the expression: 

d_ L2 = f p )(r 2 (f ■ f p ) - {f- r p )(f- f))}, Bl 
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where D is the separation distance, r is the radius vector of the star with coordinates 
(x(t), y(t), z(t)) and r p is the radius vector of the perturber with coordinates (x p (t) = 
Dcos(n),y p = -Dsin(n), 0), where n = Q p t + n and Q p = ^ tjjf. In the expression in 
the braces in (Al) we can use the unperturbed values of f(t) and r p (t). Since the orbit of the 
star is in general inclined with respect to the orbital plane of the perturber, it is convenient 
to introduce new coordinates associated with the orbital plane of the star where the radius 
vector of the star has two non-zero components: (x(t),y(t),0), and relate these coordinates 
to the coordinates (x(t),y(t), z(t)) with help of the Euler angles. We have 

x = (cos (a;) cos(x) — sin (a;) sin(x) cos(i))x — (sin (a;) cos(x) + cos(u) sin(x) cos(i))y, B2 

y = (cos(a;) sin(x) + sin(cj) cos(x) cos(i))x — (sin(a;) sin(x) — cos(u;) cos(x) cos(i))y, B3 

where i is the inclination, uj is the argument of pericentre and x is the longitude of ascending 
node. x(t) and y(t) are given by the standard Keplerian expressions 

x(t) = a(cos(£) — e), y(t) = aesin(£), B4 



where a is the semi- major axis, e is the eccentricity, e = \J(1 — e 2 ) and £ is the eccentric 
anomaly. Substituting (B2-B4) into the expression in the braces in equation (Bl) we have 

A = -eD 2 a 4 Bt Bh 

where the variable A represents the expression in the braces and 

B = (1 - ecos(0)((cos(0 - e)esin(0(a 2 - (3 2 ) + ((cos(0 - ef - e 2 sin 2 (0)«/3), B6 

a = cos(u;) cos(0) + sin(w) sin(0) cos(i), (3 = sin (a;) cos(0) — cos(cj) sin(c/>) cos(i), B7 
and = n — \- 

After averaging of (B5) over the orbital period of the star, P or b, we obtain 

A = — !— d£A = -^-(a(3)ee 2 D 2 a 4 n, B8 

Porb JO 2 



where Q = y^f - Substituting (B8) into (Bl) and taking into account the standard relation 
L = eV GMa, we obtain: 



where 



d r 15 . 9 . a .nGm „ 
-L = -- ( aP)e\-f— B9 



a(3 = ~ sin(20) cos(?) cos(2w) + j sin(2w)(sin 2 (?) + cos(20)(l + cos 2 (i))). BIO 
2 4 
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After averaging over the binary period equation (B9) must agree with equation (10). Tak- 
ing into account that the averaged value of a(3 is equal to \ sin(2o;) sin 2 (i) and using 



e = V 1 ~ £> we obtain: 



d 15 /nq/n/oX 3 ^ „ 

-e = — qeesm(2uj) sin (?) — ) Q, Bll 



dt 8 v / w y D j 

where q = m/M. It is easy to see that equation (Bll) is identical to equation (10). 

Now let us estimate the amplitude of oscillations of the angular momentum AL on time 
scale comparable to the orbital period of the perturber. Neglecting for a moment these oscil- 
lations and assuming that the eccentricity of the star due to evolution of the orbital elements 
on the 'slow' Kozai time scale is close to its maximal value, and accordingly, the angular 
momentum is close to its minimal value L m i n , it follows from the discussion above that the 
inclination i min is minimal and u ~ ±|. Therefore, we have a(5 ~ — \ cos(vm) sin(20) from 

(A10), and equation (B9) takes the form: 

d r 15 ,. ,fa\ 3 /Gm\ . . 

— L ~ -— cos(wJ I — 1 I — — 1 sm(20). B12 
Integrating (B12) we obtain: 



where 



L min + ALcos(20), 513 



AL = Y cos(vm)g (^) 2 VGMD. BU 



This paper has been produced using the Royal Astronomical Society/Blackwell Science ETgX 
style file. 
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